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Abstract 

I study the absorbing-state phase transition in the one-dimensional contact process with mobile 
disorder. In this model the dilution sites, though permanently inactive, diffuse freely, exchanging 
positions with the other sites, which host a basic contact process. Even though the disorder 
variables are not quenched, the critical behavior is affected: the critical exponents 5 and z, the ratio 
(3/v± and the moment ratio m = (p 2 ) / p 2 take values different from those of directed percolation, 
and appear to vary with the vacancy diffusion rate. While the survival probability starting from a 
single active seed follows the usual scaling, P(t) ~ t~ s , at the critical point, the mean number of 
active sites and mean-square spread grow more slowly than power laws. The critical creation rate 
increases with the vacancy density v and diverges at a value v c < 1. The scaling behavior at this 
point appears to be simpler than for smaller vacancy densities. 
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I. INTRODUCTION 



Phase transitions to an absorbing state continue to attract much attention in nonequi- 
librium statistical physics , |5j, and have stimulated efforts to characterize the 

associated universality classes. Absorbing-state transitions arise in models of epidemics, 
population dynamics and autocatalytic chemical reactions, and underly self-organized criti- 
cality in sandpile models. The simplest models exhibiting absorbing-state phase transitions 
are the contact process (CP) S] and its discrete-time version, directed percolation (DP). 

Scaling at absorbing state phase transitions has been studied extensively, both theoret- 
ically and numerically, and more recently, in experiments [3|. A central conclusion is that 
critical behavior of the DP type is generic for models exhibiting a continuous phase transition 
to an absorbing state, in the absence of any additional symmetries or conserved quantities 
. In particular, allowing particles to diffuse at a finite rate does not affect the critical 
behavior of the contact process 101 ] . 

The effect of quenched disorder on the CP, in the form of random dilution, as well as of 
randomly varying creation rates associated with each site, has been investigated in detail, 

Noest fill , and in subsequent analyses of its unusual scaling 
17, Q, 13, 20]. Little attention has been given to the effect 
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of mobile disorder, for example in the form of diffusing passive sites. An exception is the 
recent study by Evron et al. of a CP in which certain sites have a higher creation rate than 

n 

others, with the "good" and "bad" sites diffusing on the lattice [21| . These authors report 
anomalies in the off-critical behavior of the model. 

Lack of interest in the effect of mobile disorder at an absorbing-state phase transition 
presumably reflects a belief that it should not lead to any significant changes in scaling 
properties, and might in fact be irrelevant. In equilibrium, annealed disorder at fixed con- 
centration leads to "Fisher renormalization" of the critical exponents if the specific heat 
exponent a > in the pure model, and does not alter the exponents if a < 22[]. The con- 
dition for Fisher renormalization (FR), like that of Harris' criterion 23|, can be extended to 
models for which (as in the present case) a is not defined [24j , by using the equivalent rela- 
tion, dv < 2. (In the present context v corresponds to the correlation length exponent u±.) 
A subtlety nevertheless arises in extending the notion of FR to nonequilibrium systems, 
because this phenomenon requires equilibrium between the degrees of freedom associated 



with the disorder (such as vacancy positions) and the other elements comprising the system 
(e.g., spin variables). In the present case there is no energy function that could favor (via a 
Boltzmann factor e~@ H ) one disorder configuration over another; all disorder configurations 
equally probable, since the CP degrees of freedom have no influence over the vacancy dy- 
namics. We therefore speak of "mobile" rather than "annealed" disorder. (Of course, those 
disorder configurations that are more favorable to the maintenance of activity in the CP will 
survive longer and so might dominate long-time averages over surviving realizations.) 

In the FR scenario, the exponent a is replaced by —a/(l — a), whilst the other static 
critical exponents (with the exception of rj) are multiplied by a factor of 1/(1 — a). Annealed 



dynamic exponent [25|]. In the case of the contact process, under the FR scenario, the 
critical exponents (3 and u± would be renormalized by the same factor. (If we confide 
in the hyperscaling relation a = 2 — du±, we would expect that (3 — > /3/(dv± — 1), and 
similarly for the other static exponents.) It is less clear what would happen to the exponents 
governing critical dynamics, such as z or 5, or to an order-parameter moment ratio such as 



In this work I study the CP with mobile vacancies. The vacancies diffuse freely, at rate 
D, exchanging positions with nondiluted sites, be they active or inactive, which host the 
usual CP. Aside from its intrinsic interest, processes of this sort should find application in 
epidemic modelling (wandering of immune individuals in a susceptible population), and in 
population dynamics (for example, a diffusing population of a third species in a predator- 
prey system). The presence of diffusing impurities is in principle relevant to experimental 
realizations of directed percolation in fluid media p|. It transpires that the critical creation 
rate A c diverges at a certain vacancy density v c (D) < 1, so that an active state is impossible 
for v > v c . It is thus of interest to study the limiting case of A — > oo, which corresponds 
the limit of an epidemic that propagates extremely rapidly, and can be impeded only by 
fragmenting the population. 

The remainder of this paper is organized as follows. In Sec. II I define the model and 
discuss some of its basic properties as well as the results of the pair approximation. Sec. Ill 
presents the results of various simulation approaches, and in Sec. IV the implications for 
scaling behavior are discussed. 
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II. MODEL 

In the basic contact process |6| each site of a lattice is in one of two states, "active" 
or "inactive". Transitions of an active site to the inactive state occur spontaneously (i.e., 
independent of the other states of the sites) at unit rate, while an inactive site becomes 
active at rate Xn/q, where n is the number of active nearest neighbors and q the total 
number of neighbors. The configuration with all sites inactive is absorbing. The model is 
known to exhibit a continuous phase transition between the active and absorbing states at a 
critical creation rate A c ; in the one-dimensional CP A c = 3.29785(2). (Figures in parentheses 
denote uncertainties.) The transition falls in the universality class of directed percolation. 
In the CP with mobile disorder (CPMV), a fraction v of the sites are vacant, and do not 
participate in the usual CP dynamics. The vacancies hop on the lattice at rate D. In a 
hopping transition, a vacancy trades its position with its right or left neighbor. Denoting 
vacancies, and active and inactive sites by v, 1 and 0, respectively, we have the transitions 
vl — > lv and vO — > Ov at rate D/2, and similarly for vacancy hopping to the left. (From 
here on, an "inactive" site denotes one that is not a vacancy, but that happens not to be 
active.) The typical evolution shown in Fig. [T] illustrates that regions of higher than average 
vacancy density tend to be devoid of activity, and vice-versa. 

The CPMV is related to a model studied recently by Evron, Kessler and Shnerb (EKS) 



2l| . in which half the sites are "good", meaning they have a creation rate Ai, while the 
creation rate at the other, "bad" sites is A 2 < Ai. The good and bad sites diffuse on the 
lattice at rate D, independent of their state (active or inactive). The CPMV resembles the 
limit A2 — > of the EKS model, but the correspondence is not exact, since bad sites can in 
fact become active in the latter system, even if A2 = 0. (For A2 = 0, a bad site can never 
produce offspring, but it can always be activated by an active neighbor. In the CPMV by 
contrast, the vacancies can never be activated.) 

While the EKS model with A2 > represents a weaker form of disorder than the CPMV, 
one might expect the two models to fall in the same universality class. In principle, both 
models should have the same continuum limit, given by the following pair of stochastic 
partial differential equations: 



-£ = VV 2 p + (a + 7 0)p - bp 2 + 77(x, t) 



(1) 



FIG. 1: Typical space-time evolution, L = 200, v = 0.1, D = 1, and A = 4.095. Red points: active sites; 
black points: vacancies. This image corresponds to 200 time steps, with time increasing downward. 



and 

= W 2 + V-£(x,t) (2 ) 
where p(x, t) is the activity density and 0(x, t) the density of nonvacant (or "good") sites, 



rout the term oc (ftp, Eq. ([T]) is the usual 
with r](x, t) a zero-mean, Gaussian 



whose evolution is completely independent of p. Wit 
continuum description of directed percolation 
noise with autocorrelation (r/(x, t)i](y, s)) oc p(x, t)5(x — y)5(i — t'). The second term on 
the r.h.s. of Eq. ([2]) similarly represents a conserved delta-correlated noise. It is not my 
intention to analyze these equations here. They merely serve to highlight the observation 
that if, as seems plausible, they represent the continuum limit of both models, then the 
models should have a common scaling behavior. 

Suppose that, as seems reasonable, the CPMV exhibits a continuous phase transition to 
the absorbing state at a critical creation rate X c (v,D). For any nonzero vacancy concen- 
tration, A c must diverge as D — > 0, since we are then dealing with a CP on the line with 
frozen dilution, which has no active phase for finite A. The limiting value of A c as D —>■ oo 



is a subtler issue. If all three species were to diffuse infinitely rapidly, the system would be 
well mixed, and would exhibit mean-field like behavior. Note however that vacancy diffusion 
does not change the ordering of the active and inactive sites. Thus infinitely rapid vacancy 
diffusion does not correspond to a mean-field limit, but rather to a system in which a frac- 
tion v of sites are momentarily replaced, independently at each instant and at each lattice 
site, by inert elements. Hence we should expect to observe the behavior of the pure model, 
at an effective reproduction rate A' = (1 — v) A, so that the critical point is renormalized to 
Ac = X c ,pure/ (1 — v) in this limit. 

While disorder that diffuses infinitely rapidly is irrelevant to critical behavior, one may 
argue that it is relevant for finite D. Consider a correlated region in the CP, with character- 
istic size £ and duration r. If fluctuations in the vacancy density on this spatial scale relax 
on a time scale t<$, <C r, then the CP will be subject, effectively, to a disorder that is uncorre- 



lated in time, which is irrelevant 



27| . But since the fluctuations relax via diffusion, ~ £ 2 



In the neighborhood of the critical point, £ ~ |A — \ c \~ u± and r ~ £ z , with z = v\\/v±, so 
that T0 ~ r 2//z , suggesting that diffusing disorder is relevant for z < 2, provided of course 
that quenched disorder is relevant. (The two conditions are conveniently written as dv± < 2 
and dv\\ < 4.) In directed percolation this inequality is satisfied in d < 4 space dimensions. 

It is perhaps worth mentioning that the present model is not equivalent to the diffusive 
epidemic process (DEP) 28) . In the latter model, there are two kinds of particle (healthy 
and infected, say), which diffuse on the lattice. Infected particles recover spontaneously 
while a healthy particle may become infected if it occupies the same site as an infected 
particle. (There is no limit on the occupation number at a given site in the DEP.) While 
it is tempting to identify the active and inactive sites of the CPMV with (respectively) the 
infected and healthy particles of the DEP, the analogy is not valid since, in the CPMV, 
vacancy diffusion does not cause the nondiluted sites to change their relative positions, as 
noted above. 

A simple and often qualitatively reliable approach to estimating the phase diagram is 
via n-site approximations, where n = 1 corresponds to simple mean-field theory. At order 
n, the m-site joint probability distributions (for m > n) are approximated on the basis of 



the n-site distribution 
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301 ]. Since simple mean-field theory is insensitive to diffusion, 



the lowest order of interest here is n = 2 (the pair approximation), which is derived in the 
Appendix. The resulting transition line, X C (D), is compared with simulation results (for 



v = 0.1) in Fig. [2J The pair approximation turns out to be in rather poor agreement with 
simulation. In particular, it fails to show that A c diverges as D — > 0. Some improvement 
might be expected using larger clusters, but it is not clear if the correct asymptotic behavior 
would be captured even for large n. 
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FIG. 2: Critical reproduction rate A c versus diffusion rate D, for dilution v = 0.1. Points: simulation; solid 
line: pair approximation; dashed line: limiting (D — > oo) value, A c = \ c ,pure/v. 

For a fixed diffusion rate D > 0, we expect A c to increase monotonically with vacancy 
concentration v. One may ask whether A c remains finite for any v < 1, or whether it diverges 
at some vacancy concentration that is strictly less than unity. Now, for very large A, a string 
of two or more active sites is essentially immortal, for each time a site becomes inactive 
(which occurs at rate 1), it is immediately reactivated by its neighbor. Thus only isolated 
(nondiluted) sites can remain inactive. For v sufficiently large, we expect an isolated active 
site to become inactive before it can make contact with another nondiluted site, so that 
A c — > oo at some critical vacancy concentration v c (D) < 1. The following simple argument 
furnishes an estimate of v c . The rate of loss of activity due to isolated active sites becoming 
inactive is ~ pv 2 (that is, the mean- field estimate for the density of isolated active sites). On 




the other hand, the rate of reactivation of inactive sites may be estimated as Dpv(l — p — v ) , 
i.e., the vacancy diffusion rate times the mean-field estimate for the density of active- vacant- 
inactive triplets. The resulting equation of motion, dp/dt = —v 2 p + Dpv(l — p — v), has 
the stationary solution p = 1 — v — v/D, so that the activity density vanishes at v — 
D/(l + D) = v c . While this approximation is too crude to yield a quantitative prediction, 
the expression for v c does at least attain the correct limiting values for D — > and D — > oo. 
(The pair approximation does not yield a useful prediction for v c .) Simulation results for v c 
are reported in Sec. Hid. 



III. SIMULATIONS 



I performed three series of simulations of the model. The first series follows the evolution 
starting from a maximally occupied configuration (all nonvacant sites occupied) from the ini- 
tial transient through the quasi- stationary (QS) state, in which mean properties of surviving 
realizations are time- independent. These simulations will be referred to as "conventional", 
to distinguish them from the second series, which employ the QS simulation method 33]. 
The third series studies the spread of activity starting from an initial seed. In all cases, the 
initial positions of the vL vacancies are selected at random from the L sites of the ring. A 
new selection is performed at each realization. 



A. Conventional simulations 



In conventional simulations (CS), the fraction p of active sites and the moment ratio 
m = (p 2 ) I p 2 are calculated as functions of time. The time increment associated with 
a given event is At = [(1 + X)N a + DNv]" 1 , where N a and N v denote the numbers of 
active sites and of vacancies, respectively. p(t) and m(t) are calculated as averages over 
surviving realizations in time intervals that for large t represent uniform intervals of In t, a 
procedure sometimes called "logarithmic binning". In addition to the QS values of p and 
m, conventional simulations furnish the mean lifetime r from the asymptotic decay of the 
survival probability, P s ~ e~ l l T . Estimates for the critical exponents S and z are obtained 
using the relations p(t) ~ t~ 5 and m(t) — 1 ~ t x l z The latter relations describe the 

approach to the QS regime in a large system at its critical point. Finite-size scaling (FSS) 



theory implies that, at the critical point, pqs ~ L~^l Vl - , while r ~ L z , and that m converges 
to a finite limiting value as L — > oo 32J. Note that we have, in principal, two independent 
ways to estimate the exponent z, one involving the scaling of the lifetime with system size, 
the other using the approach of m to its QS value. (These scaling behaviors have been verified 
at absorbing-state transitions in various models without disorder; scaling in the presence of 
quenched disorder is generally more complicated. Whether simple FSS applies for mobile 
disorder is an open question; the results shown below provide partial confirmation.) 

I performed CS for system sizes L = 100, 200, 1600, averaging over N r independent 
realizations, where N r = 100 000 for L = 100, and decreases gradually with system size, to 
a value of 2 000 - 5 000 for L = 1600. The run times are such that all realizations eventually 
fall into the absorbing state. In analyzing the CS data, I use two criteria for determining the 
critical reproduction rate A c : (1) power-law scaling of the QS density, pqs, with system size, 
and (2) approach of the moment ratio mqs to a well defined limiting value with increasing 
system size. Verification of (1) is facilitated by studying the curvature of Id.Pqs versus 
InL, by calculating a quadratic fit to the data, and determining the range of A for which 
the quadratic term is zero to within uncertainty. It is also useful to plot L )3 ^ u± pqs versus 
L, on log scales, to check for curvature visually, as illustrated in Fig. [3j (Here (3/v±_ is 
estimated from the data for the value of A showing the least curvature.) Verification of the 
second criterion is facilitated by plotting m versus 1/L; for subcritical values this graph 
veers upward as 1/L — > and vice- versa, as shown in Fig. HI 

The behavior of the activity density as a function of time, at criticality, is shown in Fig. 
El for the case v = 0.1 and D = 1. In the pure contact process the activity follows a power 
law, with very weak corrections to scaling, for comparable times and system sizes [l|. Here 
by contrast, there is a significant positive curvature to the plots of pit) on log scales. It is 
nevertheless possible to estimate the exponent 5 from the data at longer times and larger 
system sizes. Since the graphs of p(t) are not linear (on log scales) it is not possible to 
collapse the data onto a unique master curve describing both the stationary and transient 
regimes, as is in fact possible in the contact process without dilution. The quantity m(t) — 1 
shows cleaner power-law scaling over a somewhat larger period, as shown in Fig. [6], allowing 
a rather precise estimate for the exponent z. 

Table I summarizes the results of the conventional simulations. The critical parameters 
such as /3/v± and m are calculated for each value of A simulated, and then evaluated at 
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FIG. 3: Scaled order parameter p* = L@/ U± p versus system size in conventional simulations for v = 0.1, 
D = 2, and (lower to upper) A = 3.905, 3.910, 3.915, and 3.920. Error bars are smaller than symbols. 

A c using linear interpolation. Each fit to the data involves an error bar, to which we must 
add the (typically larger) uncertainty induced by the uncertainty in A c itself. In Table I, 
the values cited for the dynamic exponent z are obtained from the growth of the moment 
ratio (using m — 1 ~ t 1 ^), and not from the FSS relation for the lifetime, r ~ L z . In fact, 
the growth of the lifetime appears to be slower than a power law at the critical point, as 
shown in Fig. [7] for v = 0.1 and D = 1. In all the cases shown (including A = 4.101, above 
the critical value), the plots of r* = L~ z t curve downward. (In this plot, for purposes of 
visualization I use z = 2.4, derived from a linear fit to the data form r(L) at A = 4.099.) 

As noted in the Introduction, averages restricted to surviving realizations may yield 
statistical properties of the vacancies that differ from those of a random mixture. I studied 
the density of nearest-neighbor vacancy pairs; for a random mixture p vv = N(N — 1)/L 2 , 
where N = v L is the number of vacancies. At the critical point (A = 4.099, D = 1, and 
v — 0.1), I find that p vv grows with time, reaching a value about 2% greater than the random 
mixture value in a system with L = 100. For system sizes L = 200 and 800, the increase 
over the random mixture value is 1% and 0.2%, respectively. Thus selection effects on the 
vacancy configuration appear to be quite weak. 





FIG. 4: Moment ratio m = (p 2 )/p 2 versus inverse system size in conventional simulations for v = 0.1, 
D = 2. and (upper to lower) A = 3.905, 3.910, 3.915, and 3.920. 

B. Quasistationary simulations 

This simulation method is designed to sample directly the quasistationary (QS) regime 



of a system with an absorbing state, as discussed in detail in 33]. The method involves 
maintaining a list of active configurations during the evolution of a given realization. When 
a transition to the absorbing state is imminent, the system is instead placed in one of the 
active configurations, chosen randomly from the list. This procedure has been shown to 
reproduce faithfully the properties obtained (at much greater computational expense) in 
conventional simulations. The theoretical basis for the QS simulation method remains valid 
in the presence of variables that continue to evolve even in the absence of activity, such as 
the vacancy positions in the CPMV, or the positions of healthy particles in the diffusive 



epidemic process 



34) 



I performed QS simulations on rings of L sites, as the conventional simulations. The 
simulation yields a histogram of the accumulated time during which the system has exactly 
1, 2,...n,..., active sites, which is used to calculate the QS order parameter pqs and moment 
ratio rags- The lifetime r is given by the mean time between attempted visits to the 




FIG. 5: Order parameter p(t) for v = 0.1, D = 1, and A = 4.099. System sizes (upper to lower) L — 200, 
400, 800, and 1600. The slope of the straight line is -0.0863. 

absorbing state, in the QS regime. Results are obtained from averages over 10 independent 
realizations, each lasting t max = 10 8 time units, with an initial portion (10 4 to 10 5 time steps, 
depending on the system size), discarded to ensure the system has attained the QS state. 



One thousand saved configurations are used to implement the QS simulation procedure [33 ]. 
Values of the replacement probability p rep range from 10~ 3 to 6 x 10 -5 (smaller values for 
larger L); during the initial relaxation period, p rep = 0.01. 

Comparing the results of QS and conventional simulations near the critical point, I find 
no significant differences for v — 0.1 and D = 0.5, in systems of up to 1600 sites. For 
D — 1, the results for p and m again agree to within uncertainty, but the QS simulations 
yield slightly smaller values for the lifetime than conventional simulations. (For L = 800, 
the difference is about 5%.) The difference can be understood if we suppose that certain 
initial vacancy configurations are more favorable than others for the survival of activity. In 
the conventional simulations each initial configuration makes a single contribution to the 
average survival time, whereas in QS simulations, realizations starting with less favorable 
configurations (with shorter lifetimes) visit the absorbing state more frequently than those 
starting from favorable ones, and so tend to dominate the average for the lifetime. This 




FIG. 6: Evolution of m(t) — 1 for v = 0.1, D = 1, and A = 4.099. System sizes as in Fig. Fhe slope of 
the straight line is 0.40. 

presumably reflects the slow relaxation of the vacancy density fluctuations. An advantage of 
conventional simulations in this respect is that many more initial configurations are sampled 
than in the QS simulations. (The former also provide information on the relaxation process, 
that is not available in QS simulations.) Despite the minor difference in lifetime estimates, 
the otherwise excellent agreement between the two simulation methods lends confidence to 
the results obtained. 



C. Spreading simulations 

In light of the surprising results of the conventional simulations, it is of interest to study 
the model using a complementary approach, which follows the spread of activity, starting 
from a localized seed. Such spreading or "time-dependent" simulations were pioneered 



in the context of contact process by Grassberger and de la Torre [35[. Starting with a 
single active site at the origin, one determines (as an average of a set of many realizations) 
the survival probability P(t), mean number of active sites n(t), and mean-square spread, 
R 2 (t) = (5^. Xj(t) 2 ) /n(t). (Here the sum is over all active sites, with Xj the position of the 
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FIG. 7: Scaled lifetime r* — L z t versus system size for v = 0.1, D = 1, and (lower to upper) A = 4.097, 
4.099, and 4.101. 

j'-th active site, and the brackets denote an average over all realizations.) The spreading 
process is followed up to a certain maximum time t m , with the system taken large enough 
that activity does not reach the boundaries up to this time. The expected scaling behaviors 
at the critical point are conventionally denoted P(t) ~ t~ s , n(t) ~ t v and R 2 {t) ~ t Za . Away 
from the critical point, these quantities show deviations from power laws. P(t), for example, 
decays exponentially for A < A c , and approaches a nonzero asymptotic value, P^, for A > A c . 

In the CP with mobile disorder, this scaling picture appears to hold for the survival 
probability but not for n(t) and R 2 (t). For example, Fig. [8] shows the decay of the survival 
probability for v — 0.1 and D = 2; P(t) is well fit by an asymptotic power law with 
5 = 0.0971(5). (Note however that at shorter times (t < 100) the decay is governed by a 
larger exponent, with a value of about 0.147.) The small differences between estimates for 8 
characterizing the decay of the order parameter in CS, and that of the survival probability in 
spreading simulations (S c and S s , respectively, in Table I), can be attributed to the deviations 
from a pure power law noted above in the CS. (Note also that the uncertainties reported 
for S s do not include any contribution due to the uncertainty in A c , and so are smaller than 
the uncertainties in 5 C .) The behaviors of n and R 2 at the critical point are illustrated in 



Fig. While I have not found a simple functional form capable of fitting these data, it is 
clear that the growth is slower than a power law. At short times n and R 2 follow apparent 
power laws, with exponents of 0.35 and 1.22, respectively, which are comparable to the DP 
spreading exponents rj ~ 0.314 and z s ~ 1.27. 

0.0 i — • — i — • — i — • — i — • — i — • — i — • — i — • — i 




2 4 6 8 10 12 14 
Int 

FIG. 8: Lower curve: survival probability P versus time for v — 0.1 and D = 2, A = 3.915, average over 
10 4 realizations. Upper curve: Scaled survival probability t s P(t) using S — 0.0971. 

D. Critical vacancy concentration 

At a fixed diffusion rate, one expects A c to be an increasing function of v, as noted 
in Sec. II; simulation results for D = 1 confirm this. (For v > 0.2 I determine A c in 
conventional simulations, using the criteria discussed above, in studies of systems with L < 
800.) To determine the critical vacancy concentration v c (D), at which A c diverges, I perform 
conventional simulations at an effectively infinite creation rate. This is done by (1) allowing 
only isolated active sites to become inactive, at a rate of unity, and (2) activating any 
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FIG. 9: Mean number n of active sites (lower curve) and mean-square distance from origin R 2 (upper curve) 
versus time; parameters as in Fig. [5] 

nondiluted site the instant it gains an active neighbor. (Thus a string of inactive sites are 
activated simultaneously when the right- or leftmost site acquires an active neighbor via 
vacancy motion.) These studies yield v c = 0.517(1) for D — 1, and v c ~ 0.4182(5) for 
D = 0.2. Fig. [10] shows \ c (v) for D = 1; the data suggest that A c ~ {y c — t>)~ L16 (see 
inset). The available results are consistent with P/v±, and the critical moment ratio m c , 
being constant along the critical line for fixed D = 1. 

I turn now to the scaling behavior at the critical vacancy density v — v c , A = oo. For 
D — 1, finite-size scaling analysis of the order parameter yields /3/v± = 0.184(20). Unlike in 
the case v = 0.1 discussed above, here the survival time follows a power-law with z = 1.98(3), 
close to the value for z obtained via analysis of the growth of m(t). The short-time scaling 
behavior p ~ t~ 5 is followed quite clearly (see Fig. [TT|) . without the strong crossover effects 
observed for v = 0.1 (Fig. [5]). Fig. [12] is a scaling plot of m — 1, using t* = t/L z , showing a 
near-perfect data collapse for three system sizes, and a clean power law, m — 1 ~ t 1/z . There 



is also good evidence for power-law scaling of the survival probability and of n and R 2 at v c , 
as shown in Fig. [13j Thus the scaling behavior at the critical vacancy density appears to be 
simpler than that observed for v = 0.1. Results for critical parameters at the critical vacancy 
density are summarized in Table II. For the two diffusion rates studied, the exponents and 
m r agree to within uncertainty. The spreading exponents satisfy the hyperscaling relation 



35[ 45 + 2rj = dz to within uncertainty. The principal inconsistency in these data is in 



regard to the scaling relation z = 2/z s . In fact 2/z s is about 8% (5%) greater than z m for 
D = 0.2 (D = 1). Studies of larger systems and a more precise determination of v c should 
help to resolve this discrepancy. 




FIG. 10: Critical creation rate A c versus vacancy concentration v for D — 1 The dashed vertical line marks 
v c . Inset: the same data plotted versus v c — u; the slope of the straight line is -1.16. 




FIG. 11: Order parameter versus time for v = 0.5176, A = oo, with D = 1. System sizes L — 796 (upper) 
and 1592 (lower). The slope of the straight line is -0.093. 

IV. DISCUSSION 

I study the one-dimensional contact process with mobile vacancies, verifying that the 
model exhibits a continuous phase transition between the active and the absorbing phases, 
as in the undiluted model. The scaling behavior is, however, considerably more complex 
than the Fisher renormalization scenario that might be expected by analogy with equilib- 
rium systems with annealed disorder. (The analogy, as noted, is somewhat flawed since the 
disorder variables are not in equilibrium with the other degrees of freedom.) Cluster ap- 
proximations appear to be incapable of providing an accurate phase diagram. Quantitative 
results on critical behavior are obtained using conventional simulations (CS) and spreading 
studies. Quasistationary simulations yield results consistent with CS, except for a slightly 
smaller survival time. 

The scaling behavior is anomalous in several aspects. First, the growth of the lifetime 




FIG. 12: Scaling plot of m — 1 versus t* — t/L z using z = 1.98. Parameters v — 0.5176, A = oo, and D = 1. 
System sizes L = 398 (+), L — 796 (x), and 1592 (circles). The slope of the straight line is 0.51. 

t at the critical point appears to be slower than a power law, as is also the case for mean 
number of active sites and mean-square spread in spreading simulations at the critical point. 
The latter finding is reminiscent of the one-dimensional CP with quenched disorder in the 
: orm of an annihilation rate randomly taking one of two values, independently at each site 
Jj], [l^]. In this case an intermediate phase arises, in which the survival probability, starting 
from a single active site, attains a nonzero value as t — > oo, but the active region grows in a 
sublinear manner with time. The spreading simulations also show a crossover from scaling 
with exponent values similar to those of DP at short times, to a different behavior at longer 
times. 

In Sec. II it was argued that A c — > oo as D —>■ 0, and that A c — > A CiPure /(l — v) in the 
opposite limit. The available data are consistent with this (see Fig. [2]); extrapolation of the 
data for v = 0.1 yields a limiting (D — > oo) value of A c = 3.65, quite near A CjPure /(l — v) ~ 
3.664. 
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FIG. 13: Spreading simulations at the critical vacancy density. Solid lines (lower to upper): P(t), n(t) and 
R 2 (t). The corresponding dotted lines show the scaled quantities P* = t s P, n* — t~ v n, and R* 2 = t~ Zs R 2 , 
using S — 0.086, rj = 0.307, and z s = 0.965. Parameters v — 0.516, A = 00, D = 1, system size L — 6000. 

For vacancy density v = 0.1, the critical exponents 5 and z, the ratio (3/vj_, and moment 
ratio m c vary systematically with the vacancy diffusion rate (Table I). With increasing D 
these parameters tend to their directed percolation values, as might be expected, since the 
limit D — > 00 corresponds to the usual contact process, with a renormalized creation rate. 
There is preliminary evidence that some critical properties are insensitive to varying the 
vacancy concentration v at a fixed diffusion rate. Preliminary results for critical behavior at 
the critical vacancy density v c suggest a simpler scaling behavior, without strong crossover 
effects, at this point. The static, transient, and spreading behavior is quite analogous to 
that observed in the basic contact process, but with a different set of critical exponents. 
This motivates the conjecture that for v > all points along the critical line A c (i>) flow, in 
the renormalization group sense, to the same fixed point. For smaller values of v one would 
then expect to observe a crossover between DP scaling and that of the new fixed point, 



associated with mobile disorder. More extensive studies are needed to verify and fill out this 
picture. 

The effects of mobile reported here are more extensive than those found by Evron et al. 
in a related model (EKS) 211 ] . These authors do not find, for example, changes in the static 
critical exponents. In comparing the two models, it should be noted that the disorder in the 
CPMV is stronger than that involved in the EKS model. The difference nevertheless raises 
an interesting puzzle regarding the universality of disorder effects. 

Intuitively, the anomalous scaling properties of the CP with mobile disorder are associated 
with fluctuations in the local, time- dependent concentration of vacancies, which relax slowly 
in the long- wavelength limit. Developing a quantitative theory of how such fluctuations 
change the critical behavior is a challenge for future work, with potential applications in 
epidemiology and the dynamics of spatially distributed populations. Investigation of the 
two-dimensional case is also relevant for applications, and of intrinsic interest for developing 
a more complete understanding of scaling. 
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Table I. Critical parameters of the CPMV with vacancy density v = 0.1. 5 C and 5 S denote, 
respectively, the exponent as determined in conventional and in spreading simulations. DP values 
from [36]. 



D 


A c P/v± 


m c 


z 




5s 


0.5 


4.375(2) 0.175(3) 


1.076(2) 


2.65(4) 


0.0745(15) 0.0765(2) 


1 


4.099(1) 0.191(3) 


1.085(2) 


2.49(1) 


0.086(2) 


0.0837(5) 


2 


3.915(1) 0.205(3) 


1.096(3) 


2.36(5) 


0.101(4) 


0.0971(5) 


5 


3.7746(10) 0.235(4) 


1.123(4) 


1.92(2) 


0.137(3) 


0.1293(5) 


DP 


3.29785(2) 0.2521(1) 1.1736(1) 1.58074(4) 0.15947(3) 


(=*e) 



Table II. Critical parameters of the CPMV at the critical vacancy density v c . 5 C and 5 S as 
in Table I. z m is the dynamic exponent as determined from the growth of m in conventional 
simulations. DP values from 1361. 



D 


v c (3/v±. m c z m 5 C 


5 S 


V z s 


0.2 
1 


0.4182(5) 0.174(6) 1.083(3) 1.95(4) 0.087(2) 
0.517(1) 0.184(20) 1.084(11) 1.98(3) 0.091(4) 


0.086(2) 
0.086(2) 


0.303(3) 0.95(1) 
0.307(1) 0.965(10) 


DP 


0.2521(1) 1.1736(1) 1.58074(4) 0.15947(3) 


(=*c) 


0.31368(4) 1.26523(3) 



Appendix: Pair approximation 



It is straightforward to construct the pair approximation to the CPMV. For convenience 
I denote vacant, inactive, and active sites by v, 0, and 1, respectively. In this approxima- 
tion the dynamical variables are the probabilities p(00) = (00), (01), etc. Note that since 
vacancies hop independently of the states (0 or 1) of the nonvacant sites, (vv) = v 2 in 
the stationary state, and during the entire evolution, if this equality holds initially. Since 
v = (vv) + (Of) + (If), we then have (If) = v — v 2 — (Of), leaving four independent pair 
probabilities. (The normalization condition, (00) + (11) + (vv ) + 2[(01) + (Of ) + (If )], permits 
to eliminate one further variable.) 

Consider, for example, the transition 00 — ■> 01. It occurs at rate A/2 provided the site 
to the right of the pair bears a particle, so the contribution to d(00)/dt due to this event is 
— (A/2)(001). (The mirror event makes an identical contribution.) In the pair approxima- 
tion, three-site probabilities are written in terms of the one- and two-site probabilities, for 
example, (001) ~ (00)(01)/(0). Proceeding in this manner we obtain 



dt K J v 



(00) 
(0) 



[A(01) + D(0v ) 



(3) 



d(01) 
dt 



(u)+ a(oom +dMM _ (01) 

2 ( J d 



1 + £_ + M 

. AW (0) 



2 V + (0) . 



(4) 



and 



d(0v) 
dt 



(lf) + 



D 



D, \ 



"(00)(0w) 

(Of) + (If) + (Of 



(0) 



(01)(lf) 



2 V ; (0) 



(5) 



d(n) 

dt 



A(01) 1 



(10) 
(0) 



D 



If) 



11) 2 + D 



(If) 
(1) 



(6) 



Integrating these equations numerically, one can determine the critical creation rate A c as a 
function of v and D. 



